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ABSTRACT 

In this paper, we discuss the use of a sensitivity equation method to compute derivatives 
for optimization based design algorithms. The problem of designing an optimal forebody 
simulator is used to motivate the algorithm and to illustrate the basic ideas. Finally, we 
indicate how an existing CFD code can be modified to compute sensitivities and a numerical 
example is presented. 
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1 Introduction 

A large number of identification, control and design problems may be formulated as infinite 
dimensional optimization problems. These problems arise in almost all fields of science 
and engineering and range in scope from inverse problems in seismology, to LQR and H°° 
control, to shape optimization in fluid/structure dynamics. See [4,7] for typical applica- 
tions. Although there are numerous approaches to solving these problems, each approach 
requires that some type of approximation be introduced at some point in the design pro- 
cess. Moreover, it is often the case that some iterative scheme is needed to solve the state 
equations (in black-box methods [3,8]) and the adjoint equations (in adjoint and “one-shot” 
methods [9]). Also, the optimization algorithm may itself be iterative. In any case, the 
development of computational methods for optimal design and control can produce several 
levels of approximation and the convergence properties of the overall algorithm are very 
much dependent on the approximations. In this paper we concentrate on the problem of 
computing accurate sensitivities for gradient based optimization algorithms. In order to 
keep the paper short and, at the same time illustrate the basic idea, we concentrate on a 
particular application. We give a brief description of the problem and use this problem to 
motivate the algorithm presented below. 


2 Optimal Design of a Forebody Simulator 

This problem is a 2D version of the problem described in [1,2,6]. The Arnold Engineering 
Development Center (AEDC) is developing a free-jet test facility for full-scale testing of 
engines in various free flight conditions. Although the test cells are large enough to house 
the jet engines, they are too small to contain the full airplane forebody and engine. Thus, 
the effect of the forward fuselage on the engine inlet flow conditions must be “simulated.” 
One approach to solving this problem is to replace the actual forebody by a smaller object, 
called a “forebody simulator” (FBS), and determine the shape of the FBS that produces 
the best flow match at the engine inlet. The 2D version of this problem is illustrated in 
Figure 1 (see [1,2, 5, 6]). 

The underlying mathematical model is based on conservation laws for mass, momentum 
and energy. For inviscid flow, we have that 
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The velocity components u and v, the pressure P, the temperature T, and the Mach number 
M are related to the conservation variables, i.e., the components of the vector Q, by 
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Figure 1. 

At the inflow boundary, we want to simulate a free-jet, so that we specify the total 
pressure P 0 , the total temperature To, and the Mach number Mo. We also set v = 0 at 
the inflow boundary. If tt/, P/, and T/ denote the inflow values of the ^-component of the 
velocity, the pressure, and the temperature, these may be recovered from To, Po and Mo by 
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The components of Q at the inflow may then be determined from (4) through the relations 
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The forebody is a solid surface, so that the normal component of the velocity vanishes, 
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utii + vn 2 = 0 on the forebody , (6) 


where m and 112 are the components of the unit normal vector to the boundary. Note that 
we impose (6) on the velocity components u and v, and not on the momentum components 
m and n. Insofar as the state is concerned, it is clear that it does not make any difference 
whether (6) is imposed on m and n or on u and v, since m = pu and n = pv and p ± 0. It 
can be shown that it does not make any difference to the sensitivities as well. 

Assume that at x = P the desired steady state flow Q = Q(y) is given as data on the 
line (called the Inlet Reference Plane) 


IRP={(z,y)\x = p,a<y<6}. 

Also, we assume here that the inflow (total) Mach number Mp can be used as a design 
(control) variable along with the shape of the forebody. Let the forebody be determined 
by the curve T - T (x),n < x < 0 and let p = ( M 0 ,T( •))• The problem can be stated as the 
following optimization problem: 

Problem FBS. Given data Q = Q(y) on the IRP, find the parameters p* = (M 0 *,r*( )) 
such that the functional 

J(p)= l -[\\Q o ,(0,y)-Q(yWdy 
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is minimized, where Qoo (z,y) = Qoo(z,y,p) is solution to the steady state Euler equation 

Clearly the statement of the problem is not complete. For example, one should carefully 
specify the set of admissible curves T( ) and questions remain about existence, uniqueness 
and integrability of “the” solution Qoo- We will not address these issues in this short note. 

Most optimization based design methods require the computation of the derivatives 
jj;Q,oo(x,y,p). These derivatives aTe called sensitivities and various schemes have been de- 
veloped to approximate the sensitivities numerically (see [3,5,10,11]). A common approach 
is to use finite differences. In particular, the steady state equation (8) is solved for p and 
again for p + Ap and then j^Qoo (z,y,p) is approximated by W'»( 1, y’P + A’)~ ( l« , ( r iy’F)i . This 
method is often costly and can introduce large errors. Another approach is to first derive 
an equation (the sensitivity equation) for Q 1 — y, p) and then numerically solve this 

equation. We shall illustrate this approach for the forebody design problem and present a 
comparison of the two methods. 

3 Sensitivities with Respect to the Inflow Mach Number 

First, we consider the design parameter Mq . Thus, we will derive equations for the sensi- 
tivity 
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The differential equation system (1) has no explicit dependence on the design param- 
eter Mqj so that equations for the components of Q' are easily determined by formally 
differentiating (1) with respect to M%. The result is the system 
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and where, through (3), the sensitivities (8) and (11) are related by 
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Note that (9) is of the same form as (1), with a different flux vector. In particular, (9) 
is in conservation form. As a result of the fact that (9) is linear in the primed variables, 
and that by (12) v', and P' are linear in the components of Q', (9) is a linear system in 

the sensitivity (7), i.e., in the components of Q'. 

Now, we need to discuss the boundary conditions for Q'. Except for the inflow condi- 
tions, all boundary conditions are independent of the design parameter M 0 2 . Thus, the latter 
may be differentiated with respect to M 2 to obtain boundary conditions for the sensitivities. 
For example, at the forebody where (6) holds, we simply would have that 

u'ni + v'ii 2 = 0 on the forebody . (13) 


Similar operations yield boundary conditions for the sensitivities along symmetry lines, 
other solid surfaces, and at the outflow boundary. Note that if instead of (6), one interprets 
the no penetration condition as one on the momentum, i.e., mni + 71712 = 0 on the forebody, 
then instead of (13) we would have that 

m'ni + n'j »2 = 0 on the forebody (14) 

which is seemingly different from (13). However, (6) and (12) can be used to show that 

m'n i + «'n 2 = p(u'n\ + v'n 2 ) + //(uni + 1 * 712 ) = p{u'n\ + 1 / 712 ) (15) 


so that, since p ± 0, (13) and (14) are identical. 

The inflow boundary conditions for the sensitivities may be determined by differentiat- 
ing (4) and (5) with respect to the design parameter Mj. Note that this parameter appears 
explicitly in the right-hand-sides of the equations in (4) and (5). Without difficulty, one 
finds from (5) that 
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4 Sensitivities with Respect to the Forebody Design Parameters 

We assume that the forebody is described in terms of a finite number of design parameters 
which we denote by P k ,k = 1 ,...,K, and that the forebody may be described by the relation 


y = *(x ] P u P 2 ,...,PK), «<*</? _ (18) 

We express the dependence of the state variable Q on the coordinates and the design 
parameters by Q = Q(<, x, y\ M 2 , Pi , P2, ■ ■ ■ , Pk)- We have already seen what equations can 
be used to determine the sensitivity of the state with respect to M 2 , i.e., for Q'. We 
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now discuss what equations can be used to determine the sensitivities with respect to the 
forebody design parameters P fc , k~ 1 , z.e., for 
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System (1) has no explicit dependence on the design parameters P*, so that equations 
for the components of Q* are easily determined by differentiating (1) with respect to P*, 
k = 1, . . K. This produces the systems, k = 1, . . . ,A, given by 
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Moreover, by (3), the sensitivities (20) and (23) are related by 
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for k = 1 , . . . , K. 

All boundary conditions except the one on the forebody also do not depend on the 
forebody design parameters P k , k = 1 For example, consider the inflow boundary 

conditions (4)-(5). Differentiating these with respect to P*, k = 1 , . . . , K yields that 


Pkl — m kl = Ufc/ = EkI = Pfc/ = Pfc/ = Uik/ = v k i = 0 , 


(25) 


at the inflow boundary. Now, consider the boundary condition (6) on the forebody. We 
have that on the forebody 
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Combining (6) and (26) we have that 
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along the forebody or, displaying the full functional dependence on the coordinates and 
design parameters, we have at a point (x,y) on the forebody, and at any time t, 

/ \ d$ 
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We can proceed to differentiate (28) with respect any of the forebody design parameters 
P k , k = The result is that, along the forebody for k = 1, . . K, 



where u, v , and their derivatives are evaluated at the forebody ( x,y = $(x)). 

If an iterative scheme is used to find a steady state solution of this system ((21), 
(25), (29)), then we assume that present guesses for the state variables u and v and their 

derivatives du/dy and dv/dy and for the design parameters Mq and P k , k = 1 K, are 

known. It follows that the right-hand-side of (29) is known as well and equation (29), the 
boundary conditions along the forebody for the sensitivities with respect to the forebody 
design parameters, is merely an inhomogeneous version of (27), the boundary condition 
along the forebody for the state. 

Let us now specialize to the type of forebodies considered by Huddleston, [5,6], t.e. 

K 

<&(x; Pi, P 2 , . . Pk) = '%2Pk<j>k(x ) , ( 30 ) 

*=1 

where «At(x), k = l, . . . , K , are prescribed functions, e.g., Bezier curves. In this case, 



Combining (29)-(32), one obtains that, at any point (x,$(x)) on the forebody and for each 
k — 1 , . . . , A' , 



For forebodies of the type (30), (33) gives the the boundary conditions along the forebody 
for the sensitivities with respect to the forebody design parameters P k , k = 1, . . . , K . It is 
now clear that, given guesses for the state variables u and v and their derivatives du/dy and 
dv/dy and for the design parameters M% and P*, k = 1,...,A', then the right-hand-side of 
(33) is known. 

5 Computing Sensitivities using an Existing Code for the State 

Suppose one has available a code to compute the state variables, i. e., to find approximate 
solutions of (1) along with boundary and initial conditions. In principle, it is an easy matter 
to amend such a code so that it can also compute sensitivities. 
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First, let us compare (1) with (9). If one wishes to amend the existing state code that 
can handle (1) so that it can treat (9) as well, one has to change the definitions of the flux 
functions from those given in (2) to those given in (10). Note that the solution for the state 
is needed in order to evaluate the flux functions of (10). 

Next, note that (9) and (22) are identical differential equations. Thus, the changes 
made to the code in order to treat (9) can also be used to treat (22). In fact, as long as 
the differentia] equation and any other part of the problem specification do not explicitly 
depend on the design parameters, the analogous relations will be the same for all the 
sensitivities. 

The only changes that vary from one sensitivity calculation to another are those that 
arise from conditions in which the design parameters appear explicitly. In our example, for 
the sensitivity with respect to , one must change the portion of the code that treats the 
inflow conditions (4)-(5) so that it can instead treat (16)-(17). The only changes needed 
to accomplish this are to the data of the inflow conditions. In the problem considered 
here, the nature (i.e. what variables are specified) of the boundary conditions at the inflow, 
and everywhere else, is not affected. Note that for the sensitivity with respect to Mq the 
boundary condition (13) on the forebody is the same as that for the state, given by (6). 

For the sensitivities with respect to the forebody design parameters, the inflow boundary 
conditions simplify to (25), i.e., they become homogeneous. The boundary condition at the 
forebody is now given by (29) or (33). Once again, the nature of the boundary conditions 
is unchanged from that for the state, and only the data that is specified is different. For the 
inflow boundary conditions, we may still specify the same conditions for the sensitivities, 
but now they would be homogeneous. The boundary conditions along the forebody change 
only in that they become inhomogeneous, (compare (27) and (33)). 

In summary, to change a code for the state so that it also handles the sensitivities, one 
must redefine the flux functions in the differential equations, and the data in the boundary 
conditions. The changes necessary in the code to account for any particular relation that 
does not explicitly involve the design parameters are independent of which sensitivity one 
is presently considering. 

The previous remarks are concerned only with the changes one must effect in a state 
code in order to handle the fact that one is discretizing a different problem when one 
considers the sensitivities. We have seen that these changes are not major in nature. 
However, there are additional changes that may be needed when one attempts to solve 
the discrete equations. In the numerical results presented below we use the finite difference 
code “PARC” (see [2,5]) to solve the state and sensitivity equations. However, thejollowing 
comments apply equally well to other CFD codes of this type. 

Since we are interested in the steady design problems, the time derivative in (1) is 
considered only to provide a means for marching to a steady state. Now, suppose that at 
any stage of a Gauss-Newton, or other iteration, we have used PARC to find an approximate 
steady state solution of (1) plus boundary conditions. In order to do this, one has to solve 
a sequence of linear algebraic, systems of the type 

(/ + A/A(Q| l n) ))Q^ +1) = (q^ + A<B(Q^)), n = 0,1,2,..., (34) 

where the sequence is terminated when one is satisfied that a steady state has been reached 
and where denotes the discrete approximation to the state Q at the time t = nAt. We 
denote this steady state solution for the approximation to the state by Q^. One problem of 
the type (34) is solved for every time step. In (34), the matrix A and vector B arise from 
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the spatial discretization of the fluxes and the boundary conditions. Both of these depend 
on the state at the previous time level. 

Having computed a steady state solution by (34), the task at hand is to now compute the 
sensitivities. We will focus on Q', the sensitivity with respect to the inflow Mach number. 
Analogous results hold for the sensitivities with respect to the forebody design parameters. 
Recall that given a state, the sensitivity equations are linear in the sensitivities. Therefore, 
if one is interested in the steady state sensitivities, instead of (9) one may directly treat its 
stationary version 

5F 'i + ^k = 0 


dx dy 


(35) 


Since (35) is linear in the components of Q', one does not need to consider marching 
algorithms in order to compute a steady sensitivity. One merely discretizes (35) and solves 
the resultant linear system, which has the form 


A'{ Q fc )Q' h =B'(Q fc ), 


(36) 


where Q' h denotes the discrete approximation to the steady sensitivity. The matrix A' and 
vector B' differ from the A and B of (34) because we have discretized different differential 
equations and boundary conditions. Note that A' and B' in (36) depend only on the steady 
state Q h and thus (36) is a linear system of algebraic equations for the discrete sensitivity Q' h . 

The cost of finding a solution of (36) is similar to that for finding the solution of (34) 
for a single value of n, i.e., for a single time step. The differences in the assembly of the 
coefficient matrices and right-hand-sides of (34) and (36) are minor. Thus, in theory at 
least, one can obtain a steady sensitivity in the same computer time it takes to perform one time step 
in a slate calculation. If one wants to obtain all the sensitivities, e.g., K + 1 in our example, 
one can do so at a cost similar to, e.g., K + 1 time steps of the state calculation. This 
is very cheap compared to the multiple state calculations necessary in order to compute 
sensitivities through the use of difference quotients. 

In practice, these “optimal” estimates of speed up are rarely achieved. Moreover, it is 
important to note that finite difference (FD) and sensitivity equation (SE) methods do not 
necessarily produce the same results. Since the ultimate goal is to find useful and cheap 
gradients for optimization, the most important issue is whether or not the SE method 
combined with an optimization algorithm produces a convergent optimal design as fast as 
possible. We have tested this scheme on the forebody design problem with excellent results. 


6 Numerical Example 

In order to illustrate the use of the SE method in computing sensitivities, we used the 
PARC code as described above to find approximate solutions of the sensitivity equations 
and compared the results to the finite difference method. In Figure 2 we show the ap- 
proximations of y, Mq , Pi, P?) for Mo = 2, Pi = .1 and P 2 = -15 where the forebody 

is described by two Bezier parameters (Pi,P 2 )- Both pictures are “converged” estimates. 
Note that there are considerable differences between the FD method and the SE method. 
Moreover, in Figure 3 we see that not only do the FD and SE methods produce different 
sensitivities, the value of the step size APi can greatly influence the FD approximations. 
Finally, we note that the SE method ran 4 to 5 times faster than the FD method. Also, 
although space prohibits a discussion of the optimization problem here, we have used the 
SE method in a trust region optimization scheme to produce an optimal forebody design 
for the 2D problem in [5,6]. These results will appear in a forthcoming paper. 
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7 Conclusions 

The problem of computing accurate sensitivities in problems involving solutions to param- 
eterized partial differential equations is an important part of optimal design. The goal is to 
find derivatives of solutions of partial differential equations with respect to various param- 
eters (including domain shapes) and to use these derivatives in some type of optimization 
scheme. In almost all practical problems, solutions must be obtained by numerical approx- 
imations. This fact leads to “black box” methods for optimal design. In its most basic 
form, a black box method produces approximate solutions that are then differentiated (by 
finite differences, automatic differentiation, etc.). The sensitivity equation (SE) method 
presented here is based on first deriving partial differential equations for the derivatives 
and then approximating these equations numerically. Both approaches produce numerical 
approximations of the sensitivities. However, the (SE) method can often reduce computa- 
tional effort, speed up the calculations and, provided that accurate computational schemes 
can be devised for the sensitivity equations, the derivatives can be computed with the same 
degree of accuracy as the state. The 2D optimal forebody simulator problem is an excellent 
problem for illustrating these points. The numerical results presented here show that the 
(SE) method is potentially applicable to real problems and, at the same time, raises many 
interesting theoretical and practical questions. 
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Sensitivity of U-Momentum 
with respect to First Bezier Parameter 


Sens itivity Equat ion Me thod 



Finite Difference Method 



Figure 2 


Values taken at: 
Ini. Mach # = 1.7 
Bez. P. #1 = 0.10 
Bez. P. #2 = 0.15 


Legend 
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Absolute Difference of U-Mom. Sensitivities 
with respect to First Bezier Parameter 
obtained using F.D. and S.E. Methods 


Values at: Inlet Mach # = 1.7 

Bezier Parameter #1 = 0.10 
Bezier Parameter #2 = 0.15 



Step Size=0.0001 Step Size=0. 00001 


Figure 3 
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